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Abstract — This paper proposes the use of quadrangular Bezier 
patches to model and modify the shape of the boundary for 
linear elasticity problems in 3D. The representation of the 
boundary in this way derives directly from computer graphics 
and have been analytically included in developed by the 
authors parametric integral equation systems (PIES). PIES 
are the modified classical boundary integral equations (BIE) 
in which the shape of the boundary can be modeled using a 
wide range of parametric curves and surfaces. The proposed 
approach eliminates the need for domain and boundary 
discretization in the process of solving boundary value 
problems, in contrast to popular traditional methods like FEM 
and BEM. On the basis of the proposed approach, a computer 
code has been written and examined through numerical 
examples. 

Index Terms — parametric integral equation systems (PIES), 
boundary integral equations (BIE), Navier-Lame equations, 
linear elastics in 3D, quadrangular Bezier patches 

I. Introduction 

Finite (FEM) and boundary (BEM) element methods in 
the process of solving boundary value problems replace 
continuous domain of governing equations by discrete mesh 
of finite [1,2,3,4] or boundary elements [5,6,7]. In order to 
improve the accuracy of results, it is usually necessary to 
break the domain into a larger number of smaller elements or 
use higher order elements. Additionally, numerical stability 
of solutions depends significantly on the adopted 
discretization scheme. A major focus of our research is the 
development a new method for solving boundary value 
problems based on proposed parametric integral equations 
systems (PIES), which in contrast to traditional methods 
(FEM, BEM) do not require domain or boundary discretization. 
The shape of the boundary in PIES is mathematically defined 
in the continuous way by means known from computer 
graphics parametric curves in the case 2D and parametric 
surfaces in the case of 3D boundary value problems, and its 
definition is practically reduced to giving a small set of control 
points. So far, however, PIES has been mainly used to solve 
2D potential boundary value problems modeled by partial 
differential equations such as: Laplace [8], Helmholtz [9], 



Poisson [10], and Navier-Lame [11]. High accuracy and 
effectiveness of PIES for those problems has been 
encouraging its generalization to 3D boundary problems. In 
the case of 3D problems, we have restricted our 
considerations only for Laplace [12] and Helmholtz [13] 
equations. We have also begun preliminary work 
investigating the use this strategy for the Navier-Lame 
equations to solve 3D linear elastic problem in polygonal 
domains [ 14] . The aim of this paper is to explore the potential 
benefits of using quadrangular Bezier surface patches to 
model and modify the curvilinear shape of the boundary in 
3D problem of linear elasticity modeled by Navier-Lame 
equations and solved by PIES. The formula of PIES for the 
Navier-Lame equations has already been obtained by 
adopting the general assumption that the boundary can be 
described in general way by any parametric functions. 
However, the assessment of efficiency of using Bezier patches 
for boundary representation requires direct substitution them 
into PIES. This has been done on the basis on written 
computer program whose results are presented in the last 
part of this article. 

II. QUDRANGULAR BEZIER PATCHES 

In obtained PIES for 2D boundary problems, the boundary 
geometry is considered in its mathematical formalism and 
defined with help of parametric linear segments [11] and also 
by parametric curves [8,9]: Bezier, B-spline and Hermite. The 
proposed strategy for 2D boundary description has been 
proved to be effective both from the point of view of 
simplifying the boundary declaration as well as improving 
the accuracy of obtained solutions compared with FEM and 
BEM. Boundary problems in 3D are characterized by much 
higher complexity. In this article we evaluate quadrangular 
Bezier surface patches as boundary representation in 3D. 
Quadrangular Bezier patches are an extension of parametric 
Bezier curves. Both curves and surfaces can be simply 
declared by a limited set of control points. Quadrangular Bezier 

surfaces are specified by an array of n x m control points Py 
and written mathematically by [ 1 5] [ 1 6] 
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p(v, W ) = ££^B fm 

i=0 j=0 



< v, w < 1, 



(1) 



where Z? im (v) , (w) are Bernstein basis functions written 
in the form analogous like for Bezier curves 



Bj v )= m v^l-vr bJw)- 



'(1-w)"- 



(2) 



The main advantage of Bezier patches lies in the simplicity of 
their creation and modification with help of a limited number 
of control points. In this paper, we deal with Bezier patches 
for n = m = 1 and n = m = 3 (also referred here as degree of 
1st and 3rd). Bezier patch of first degree takes the form of a 
quadrangular plate defined by four points ^oo,^o>^n>^oi 
placed at the corners of defined rectangle as shown in Fig. 1 . 

1. - 




Figure 1. Quadrangular Bezier surface patch of degree 1 with 4 
corner points. 

The patch is described by following matrix form [17 



P(v,w) = [l-v v] 



^00 




1 — w 




V 


w 



< v, w < 1 



(3) 



In turn, Bezier patch of degree 3 is defined by 16 control 
points. Its visualization is shown below. 




J5. 



Figure 2. Quadrangular Bezier surface patch of degree 3 with 16 
control points 

The patch of degree 3 can be described by following matrix 
form [17] 



P(v, 



M ■■ 
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W 








y. 



(4) 



Fig. 3 show sample modification of Bezier patch of the third 
degree by moving individual control points. 



b) 



■ 



* 



■Vji 



Figure 3. Examples of shapes defined with Bezier patch of degree 3: 
a) the initial flat surface, the updated surfaces after moving selected 



III. 



b) Poo > c ) PlO an d d) P()3 control points. 

Modeling And Modification Of The Curvilinear 
Boundary Surface By Control Points 
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A single Bezier patch can only capture a small class of 
shapes. The potential ability to model more complicated 
boundary shapes lies in the increase the degree of the 
polynomials used in the basis functions (2). From practical 
point of view it is connected with increase the number of 
used control points. However, the relative lack of local shape 
control by patches of higher degree makes them more difficult 
to accurate capture the real shapes. In practice, the most 
common solution is to join independent low degree patches 
together into a piecewise surface. Fig. 4 shows declaration of 
the boundary in PIES defined by five Bezier patches of first 
and the one of the third degree. The boundary is practically 
defined by 20 control points, of which 16 belong to curvilinear 
patch of degree 3. By appropriately selecting and moving 
one or more control points it becomes possible to effective 
change the shape of the declared boundary. In many inverse 
problems, concerning for example shape identification or 
optimization, it is necessary to perform multiply modifications 
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of initial shape of the boundary. The declaration of boundary 
geometry on the basis of Bezier patches allows for easy 
modification of the boundary by moving individual or group 
of control points. Fig 4b shows the modification of the initial 

domain by moving only one P 22 point. 



a) b) 




Figure. 4. Modeling the boundary by Bezier patches: a) posed 
control points, b) modification of initial shape by moving P 2 2 ■ 

Two variants of the geometry from Fig. 4 show the 
effectiveness and simplicity of proposed boundary shape 
modification performed using a small input data by only 
control points. Fig. 5, in turn, shows the further steps to 
create a part of a hollow cylinder with flat and curvilinear 
Bezier patches of degree 1 and 3 respectively. 

3) V) 




Figure. 5. Further stages of defining a part of a hollow cylinder 
bounded by Bezier patches: a) the patch of the first degree, b), c) 
patches of the third degree, d) the complete geometry. 

The domain is described unequivocally and accurately only 
by 48 control points of the individual surface patches. The 
calculation required to find the locations of control points 
for Bezier patches of degree 3 from the example is presented 
in Fig. 6. 
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Figure. 6. Determination of the exact locations of the control 
points for Bezier patch from the base of the cylindrical domain. 

The outer edge of the patch from Fig. 5 formed by the points 
^oo ' ^oi ' ^02 ' ^03 snou ld take the shape of a quarter circle of 
radius r ■ Assuming that the center point of the circle is 
(0,0,0) and radius r — \ the coordinates of control points 

are as follows P 00 (1,0,0) , P m (l,k,0) , P Q2 (kX0) , 
P 03 (0,1,0) . The coefficient k may be determined by formula 
(5), where for v - , w = 0.5 the outer edge of the patch 

should pass through the point ( V2 / 2, V2 / 2,0) , hence 

4l / 2 = (l - 0.5) 3 + 3*(l - 0.5) 2 0.5 + 3(l - 0.5)0.5 2 + 0.5 3 (5) 
and it = 4(V2-l)/ 3. 

The location of other control points from boundary geometry 
from Fig. 5 are determined in the same way, taking into account 
the appropriate value of the radius. Presented boundary 
representation by Bezier patches is directly included in 
presented in the next section mathematical formalism of PIES. 
Therefore, modeling the boundary is only limited to the above- 
mentioned declaration of surface patches without their further 
division, for example into boundary or finite elements. This is 
a significant advantage of PIES, which leads to simplification 
of both the boundary geometry description and numerical 
calculations. 

IV. PIES For The Navier-Lame Equation 

In the previous section, we have shown that the shape of 
3D domains can be modeled using parametric surface patches. 
This is very effective from the point of view of interactive 
modification and creation of 3D domains and also their 
visualization. Although in the boundary value problems we 
also describe the shape the domain and boundary, the main 
focus of our attention is not the visual representation of their 
geometry, but finding solutions on the boundary or in the 
domain for given boundary conditions, described shape of 
the domain and the introduced differential equations that are 
solved. This makes the problem much more complex. In 
boundary value problems solved using traditional numerical 
methods it is needed to divide computational continuous 
physical domain into elements. This approach requires a large 
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amount of elements, but it is easy to apply to the classical 
boundary integral equations (BIE). BIE is characterized by 
the fact that the boundary is defined in general way by the 
boundary integral [18]. Therefore, BEM used to solve BIE 
requires discretization of the boundary. Such boundary 
representation has also its disadvantages, because any 
modification of the boundary shape requires a re- 
discretization of the boundary. Unfortunately, regardless of 
its efficiency presented in the previous section the way of 
modeling and modification of the curvilinear boundary surface 
by parametric patches is not possible to apply directly in 
classical BIE. Therefore an analytical modification of BIE have 
been performed. This modification has been performed on 
the basis of assumptions that the boundary is not defined 
using the boundary integral but with help of discussed in 
previous sections surface patches. PIES for the Navier-Lame 
equations in 3D has been obtained as a result of the analytical 
modification of the traditional BIE. Methodology of this 
modification for 2D problems was presented in [11]. 
Generalizing this modification for Navier-Lame equations in 
3D, we obtain the following formula of PIES [14] 



v-i f f — * 

0.5m / (v 1 ,w 1 ) = 2 j I I {£/j/(v 1 ,w 1 ,v,w)/> / (v,h') 

J=1 "j-i W H 

- P t * (vj , w 1 , v, w)Uj (v, w) }Jj (v, w)dvdw, 



(6) 



where 



v 



<V<Vj 



Wj_ 1 <w<w jt i = 1,2,3..., n, and n is the number of 

parametric patches that create the domain boundary in 3D. 
These patches is included in the kernels (integrand functions) 
of the obtained PIES. Therefore PIES is not defined on the 
boundary as BIE, but on the mapping plane. Such a definition 
of PIES has the advantage that the change of the boundary 
is automatically adjust of PIES to the modified shape without 
any additional re-discretization as in the case of classical 
BEM. In this paper the patches takes form of Bezier surfaces 
that create closed 3D domain. 

The integrands U y* , Py* in equation (6) are represented in 
the following matrix form 



U ij (v[ , xv l , v, w) = 



Py*(v l ,w l ,v, xv) = 



(1 + v) 



2,7t{l-v)Eri 



1 





U n 


U u 


u 2i 


Un 


U 23 


U 3i 


u 32 


u 32 


Pn 


Pn 


Pn 


P 21 


P22 


P23 


P31 


P32 


P33 



(7) 



(8) 



The individual elements in the matrix (7) can be expressed in 
an explicit form as 



If Un =W U t -™ 



t/ 11= (3-4v)+ , 
V 
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'12 2 



'13 2 

V 



u 



V2V1 



21 " 2 

V 



U =(3-4v) + 



u 



23 



^3 



</31=» t/ 32 =^ ^33=(3-4v) + % 

while in the matrix (8) by 

Ai=((l-2v) + 3 4)^, 
77 on 

Pi2=3 m^_ (1 _ 2v) ¥2-w 



, 2 dn 



, |3=3 M^_ (1 _ 2v) »-w 



j 2 dn 



P2l=3 a^- ( i-2v)W-^ 



7 2 dn 



ri2^Sri 



P 2 2=((l-2v) + 3^) 



n 



2 ' Dn ' 



1 2 dn 



P3l=3 a^. ( i-2v)»-w 



7 2 dn 



Pi2=3 ^jL. {l . 2v) n^-^n3 



, 2 dn 



P 33 =((l-2v) + 3^)- 



4>— 

?7 2 dn 

Created by Bezier patches graphical representation of the 3D 
domain is directly incorporated into the mathematical formula 
of PIES without requiring an additional its division, for example 
into boundary elements. Kernels (7,8) include in its 
mathematical formalism the shape of a closed boundary, 
created by means of appropriate relationships between Bezier 

patches (I, j = 1,2,3..., n) , which are defined in Cartesian 

coordinates using the following relations 

7?! = Pf (V, W) - P { m (V! , W, ) , 772 = ^ (v, W> - /f > (Vj , W, ) 
(4) 

773 = Pf ( v, w) - P/ 3) (y 1 ,w 1 )md r , = [tf + n \ + n \ f 5 (9) 

where Pf\v,w), Pj 2> (v,w), Pf\v,w) are the scalar 
components of Bezier patches, described by formula (1). This 
notation is also valid for the patch labeled by index / with 

parameters Vj w x , ie. for j = I and for parameters v = v x and 

w = . It is possible to insert into (9) Bezier patches of any 
degree expressed by (1). In the paper the patches of degree 
1 and 3 is used. 
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The smooth character of Bezier surfaces requires new ways 
about computing the Jacobian J j\y,w) and normal 

derivatives n,(v, w), n 2 (v,w), « 3 (v, w) in (7,8). In the case of 
Bezier patch of degree 1 the Jacobian and normal derivatives 
are constant on the whole surface area regardless of the 
actual values for parameters v, w . While in the case of the 
patch of degree 3 these parameters change at any point of 
the patch, depending on changes of v, w . A formula for the 
Jacobian and normal derivatives can be determined 
analytically from expression ( 1 ) that define quadrangular Bezier 
surfaces as follows 



n m (v,w) 


A (v. w ) 
= m ) ' >, m = 1,2,3 
Jj(y,w) 


(10) 




where 








Ai(y..w] 




dP^iv, 


w) ap^(v,w) 


dw cH' 


Sv 


dw 




6p/ 3> (r : w)aP ; a >(v iW ) 




w) 8P^(v,w) 


dw dv 


dv 


dw 






6P/»(v 


w)eif 2) (v,w) 


dw dv 


dv 


■dw 


and 









J j (v, w)= [A* [v, w)+ A\ (v, w}+ A* (v, w)f 5 . 
The outward or inward direction of normal derivatives 
H[(v, w), n 2 {v,w), n 3 (v, w) is determined by the appropriate 
numeration of control points of Bezier patches. The use of 
PIES for solving 2D and 3D boundary problems made it 
possible to eliminate the need for discretization of both the 
mentioned above boundary geometry, and the boundary 
functions. In previous research works, the boundary 
functions both defined as the boundary conditions, as well 
as obtained after solving PIES are approximated by Chebyshev 
polynomials [8]. The numerical solution of PIES for 3D 
problems is very similar to that previously described for 2D 
problems [9]. The difference appears in the use of more 
complex series to approximate the solutions, which are 
presented in the following form 

N M 

u j ( V ,w) = YZ u T )T J P) ^ r j r)( ^, (11) 

p=Q r=0 
N M 

pj "0 = Z Z pj pr)T i p) w) r \ W ) , (12) 

p=0r=0 

where Uj ,pF are requested coefficients and rj p) (v), 
Tj r) (w) are Chebyshev polynomials. One of these functions, 

either Uj(v,w) or Pj(v,w) , depending on the type of the 

resolved boundary problem, is posed in the form of boundary 
conditions, whereas the other is the searched function 
resulting from the solution of PIES. Both the posed and 
searched boundary functions are approximated using the 
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generalized approximating series (11) or (12). Finally, after 
substituting these series to formula (6) and writing it at the 
collocation points, we obtain an algebraic equation system 

with respect to the unknown coefficients u P or Pj . 

V. Numerical Examples 

In order to evaluate proposed PIES for defined by Bezier 
patches domains shown in Fig. 4a and 5d, a numerical 
verification was performed with displacement Dirichlet 
boundary conditions for the Navier-Lame equations obtained 
from equations (13). These formulas are also exact solutions 
of the Navier-Lame equations [19] 

i h (x) = 1 ^ 3 , u 2 (x)= 1 =g 3 , (13) 



The normal directives of functions (13) computed on the 
surface of the Bezier patches can be identified with exact 
analytical solutions of the problem on the boundary. The 
derived normal derivatives for these functions are presented 
below 

Pi (x) = = 2 ni (x) + 0An 2 (x) + 0.4n 3 (x) , 

dn 

p 2 (x) = iM*) = o.4n, (x) + 2« 2 (x) + 0.4« 3 (x) , (14) 
dn 

p Jx)= _ QArh ( X ) + o.4« 2 ( X ) + 2n 3 (x) . 

dn 

Table 1 shows the accuracy of solutions obtained on the 
boundary formed by Bezier patch of degree 3 for the domain 
shown in Fig. 4. We have analyzed different shapes of this 
patch, modified by changing the position of four central 

control points P n , P l2 , P 2l , P 22 . We choose E = \Pa and 

v=0.25 for all calculations. 



Table. 1 . The errors norm L 2 of solutions in PIES for various shapes 
of Bezier patch of degree 3 from Fig. 4. 



C ontrol points 


Ths errors norm 




I^lj,^ 


-Fh (0.33, 0.33, 1.25) 




fli (0.33,0.67, 1.25) 


0.0492533 % 


0.0507693 % 


ill (0.67,0.33, 1.25) 


0.0959662 % 


1*22, (0.67,0.67,1.25) 




ill (0.33, 0.33, 1.5) 




P_ 2 ::::\:o.::\r. 


0.100674% 


0.0S6942 % 


ill (0.^7,0.33,1.5) 


0.119345 % 


il2 (0.67,0.67,1.5) 




ill (0.33, 0.33, 1.5) 




H 2 (0.33,0.67,1.5) 


0.193235 % 


0.260095 % 


ill (0.67,0.33,1.5) 


0.152921 % 


P 22 (0.67,0.67,0.5) 
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The error norms L 2 for solutions from column 2 are calculated 
on the basis of the following formula [20] 



'1 



,(*) 



= 100%, 



(15) 



where p^' represents a set of K = 85 solutions in PIES 
obtained on the boundary represented by the patch placed 

on upper face of the cube, while are exact solutions 

given by equations (14). The solutions on the boundary were 
obtained with N = M = 9 coefficients in the approximating 
series (12). Table 2, in turn, lists the errors obtained in PIES 
related to exact values ( 14) at selected points on Bezier patch 
of degree 3 which forms an outer cylindrical area in Fig. 5d. 

Table. 2. The relative error of solutions at selected points on outer 
Bezier patch of degree 3 from Fig. 5d. 





Relative error 


Xj,X 2 ,X3 


Pl,P2,P3 [%] 


5.32562. 1.44137.1.39236 


0.561040.0335302.0.559113 


3.16054. 5.10106. 139236 


4303 1 7 . 019696. 07 1019 


5.52349. 2.33547. 1.39236 


0.103074. 0.429155. 0.262713 


2.33547, 5.52349.0.39616 


0.304394,0055514,0.125395 


1.441S7.5.S2562, 0.99251 


0.343432. 0.000305. 0.SS944S 



The error norms L\ (15) for solutions obtained in 340 points 
on all six Bezier patches from geometry shown Fig. 5d for 
components Py,p 2 ,Pj, take values 0.40159%, 0.415836% and 
0.700875 % respectively. Low level of this errors demonstrate 
a sufficiently good agreement of solutions in PIES with exact 
solutions (14). 

Conclusions 

In this paper, our purpose is to discuss the possibility of 
using quadrangular Bezier surface patches with different 
degrees to describe 3D geometry for solving linear elastic 
problems in PIES. This approach in conjunction with PIES 
makes it possible to develop a completely non -element method 
for solving 3D boundary problems. Numerical example 
demonstrates the simplicity of the modeling procedure and 
the good accuracy of the results. 
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